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Using the fractional interaction law to model the impact dynamics in arbitrary form 

of multiparticle collisions 

Jacek S. Leszczynskf] 

Czestochowa University of Technology, Institute of Mathematics and Computer Science, 
ul. Dabrowskiego 73, ^2-200 Czestochowa, Poland 

Using the molecular dynamics method, we examine a discrete deterministic model for the motion 
of spherical particles in three-dimensional space. The model takes into account multiparticle colli- 
sions in arbitrary forms. Using fractional calculus we proposed an expression for the repulsive force, 
which is the so called fractional interaction law. We then illustrate and discuss how to control (cor- 
relate) the energy dissipation and the collisional time for an individual particle within multiparticle 
collisions. In the multiparticle collisions we included the friction mechanism needed for the transi- 
tion from coupled torsion-sliding friction through rolling friction to static friction. Analysing simple 
simulations we found that in the strong repulsive state binary collisions dominate. However, within 
multiparticle collisions weak repulsion is observed to be much stronger. The presented numerical 
results can be used to realistically model the impact dynamics of an individual particle in a group 
of colliding particles. 

PACS numbers: 45.05.+X, 45.70.-n, 45.50.Tn, 83.10.Pp, 83.10.Rs, 45.10.Hj 



I. INTRODUCTION 

Nature of particulate flows offers the physics and en- 
gineering communities an opportunity to analyse the in- 
teresting behaviour of granular materials. From a phe- 
nomenological point of view such a flow, being halfway 
between a solid and a liquid state, is not well understood 
because the basic physics is extremely complex. In a lo- 
cal state, the simplest form of the granular dynamics is 
as follows - during an arbitrary extortion particles move 
individually and during collisions particles may exchange 
their momenta and energies. Therefore the collision pro- 
cess plays a dominant role in the development of theoret- 
ical studies and also in the performance of simulations. 
For an understanding of the collision process we need to 
consider a simple situation, focusing on what happens 
when two particles collide. In other words, we need to be 
able to distinguish the following basic phenomena: static 
contact @ , cohesion [p f3lj| , attrition , erosion 
and fragmentation |14|. These phenomena may occur 
simultaneously or respectively when an individual parti- 
cle impacts with another. After impact separation p^j 
or clusterisation of the two particles occurs. In ad- 
dition, the particles may gain or loss mass. Here we 
will focus on the dynamics of the collision process which 
may be decomposed into impact and contact processes. 
However, as the contact process is formed, we can also 
notice rebound |24j or static contact @, or permanent 
contact, called cohesion @- These processes exist simul- 
taneously when we analyse the dynamics of colliding par- 
ticles. With regard to the granular dynamics involving 
many particles in motion, we can observe multiparticle 
collisions |29|. especially when particle concentration is 
very dense, because collisional times between several bi- 
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nary particle contacts are higher in comparison to their 
separation times. Multiparticle collisions occur when an 
individual particle collides with neighbouring particles, 
so that those contacts have direct a influence on each 
other. Only, an infinitcsimally short collisional time is 
required for binary collisions |24j. In all the considered 
cases the collision process between the two particles is 
characterised through the collisional time, which is de- 
pendent on the impact energy and the physical properties 
of the contacting surfaces. Moreover, after impact dissi- 
pation of energy occurs between the colliding particles. 
Therefore the simulations of such dynamics are limited 
by assumptions concerning the collision process. One of 
the major aspects which needs to be taken into account 
in the simulations is how to control (correlate) the colli- 
sional time and the energy dissipation associated with an 
individual particle during the dynamics of multiparticle 
collisions. 

Generally two different ways exist to model the dynam- 
ics of a granular material. The continuum approach [j| is 
based on binary collisions of smooth spherical particles. 
Unfortunately, the introduction of real quantities such as 
distribution of particle dimensions, particle shapes, their 
surface wetness and roughness, etc., greatly limit the ap- 
plication of continuum models. Balzer et al Q inform 
us that the kinetic theory is useful for the modelling of 
gas-solid flow applications in industry: where the geom- 
etry involved is complex (many different inlets or/and 
outlets). However, the kinetic theory cannot reflect the 
real dynamics involved in multiparticle collisions because 
the collisional time is defined only for binary collisions. 

The discrete deterministic approach more realistically 
reflects the collision process. Note that multiparticle 
collisions in the discrete approach are decomposed into 
several binary collisions. In this approach one may dis- 
tinguish two general methods. The molecular dynam- 
ics method |27j takes into account an expression for the 
repulsive force acting between a pair of contacting par- 
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tides. In this method particles virtually overlap when 
a contact occurs. The overlap reflects the quantitative 
deformations of particle surfaces because the modelling 
of realistic deformations would be much too complicated. 
The interaction laws [J, lla . 1 3 a] in the molecular dynam- 
ics method define basic models of the repulsive force for 
two colliding particles. They are valid for particle colli- 
sions which are independent from one another. The next 
method, called the event driven method 0, assumes 
instantaneous changes in the direction and value of par- 
ticle velocities according to conservation equations each 
time a binary contact occurs. As shown in [17| the basic 
difference between the event driven and molecular dy- 
namics methods is the collisional time between a pair of 
particles. In the event driven method this time is ide- 
ally zero. Note that this situation is quite different for 
the molecular dynamics method, where the contact time 
is greater than zero and is dependent on parameters de- 
scribing the structure of two contacting surfaces, and is 
of course dependent on the impact energy. However, the 
repulsive force models in the molecular dynamics method 
underestimate the energy dissipation in multiparticle col- 
lisions [l8l . l25| (This is the so called "detachment effect" ) 
but in the event driven method an inelastic collapse |22( 
occurs. 

In this paper we will focus on the molecular dynamics 
method because this gives us a chance to correlate the 
collisional time and the energy dissipation during multi- 
particle collisions. We shall introduce a novel mathemat- 
ical description of this method taking into account the di- 
vision of the collision process into an impact phase, a con- 
tact phase and another phase occurring after the contact 
phase. We assume that the impact phase and the phase 
formed after the contact phase are infinitesimally short 
in time. Consequently, we will analyse the well-known in- 
teraction laws of the repulsive force in the contact phase 
in order to examine several difficulties within the colli- 
sions. On the base of preliminary results [15j we shall 
introduce a novel form of the repulsive force defined un- 
der fractional calculus |2^| . We will also demonstrate the 
basic properties of this force and focus on what happens 
with the collisional time and the energy dissipation for 
multiparticle collisions. This analysis is necessary in com- 
putational simulations of the cluster dynamics. Within 
the cluster one may notice non-permanent contact and/or 
cohesion phenomena between several pairs of colliding 
particles. 



II. THE DISCRETE MODEL OF MOTION OF 
FOR INDIVIDUAL PARTICLE 

Let us turn our attention to a set of spherical particles 
moving under arbitrary extortion. The spherical shape 
of the particle makes only the mathematical description 
easier and does not make the model in any way poorer. 
The reader may find in ^20] more information concerning 
the molecular dynamics technique adapted to arbitrary 



form particle shapes. The particles are numbered by the 
discrete index i — 1, . . . , np, where np is the total num- 
ber of considered particles. We describe an individual 
particle through its radius (or diameter di), mass m;, 
inertia moment Ji , position x^ of the mass centre, linear 
speed Xi and angular velocity uji. With regard to the 
collision of two individual particles we also introduce the 
natural function ^ i by assumption) of a par- 

ticle i in order to find the particle index of a particle in 
a set of particles np. Several papers 0, 0, ^ present 
different algorithms that detect particle collisions, being 
dependent on their shapes, and consequently that to find 
the natural function For a binary collision we ne- 

glect phenomena which cause a change in the mass of 
an individual particle. Thus in our discrete model we 
do not take into account fragmentation, attrition and 
erosion which eventually take place during the collision 
process. These phenomena will be the subject of future 
investigations. 

However, after the contact, which is the second phase 
of the collision process, rebound, non-permanent contact 
(static contact) or cohesion can arise simultaneously. In 
this paper, we will try to model above the phenomena by 
introducing a novel mathematical description and a novel 
form of the repulsive force into the molecular dynamics 
method. 



A. Mapping local coordinates onto global ones and 
denning the overlap 

Starting from the description shown in Fig.^ let us in- 
troduce several definitions before formulating the motion 
equations. First, we assign local coordinates as (£, rj, () 
and global ones as (x, y, z). When we consider a contact 
which eventually takes place between two particles then 
the normal unit vector e^^) that connects the particle's 
centres of mass reads 



£ i(0 



e O'«> e Ci«' e O-(i) 



(1) 



where 



represents a norm calculated from the rela- 



tive coordinate XjU) — x,;. Tangential unit vectors which 
operate on a tangent plane (rotated by 5 to the normal) 
become 



.0 



\x,y 
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(3) 



where 



represents the norm which is calculated 



only in the tangent plane. When a particle hits a wall 
we redefine unit vectors J5J and 10 putting xb n (^ 
instead of ^-jU) , where xb^) is a point whose coordinates 
issue from the line that crosses the particle's centre of 



3 



mass and is perpendicular to the wall. The general form 
of the base vectors is presented as follows 



(4) 




FIG. 1: Scheme illustrates particle collisions with useful no- 
tations. 



Fig. ^presents the "virtual overlap" CjU) of two par- 
ticles experiencing a contact. With regard to £j we define 
the overlap as 



(5) 



which is associated with the particles having spherical 
forms. Note that only positive values of formula JSJ indi- 
cate a contact while negative ones confirm that the con- 
sidered particles are in separation. This means that they 
move individually. As presented in previous section and 
in Fig. ^ the overlap reflects the penetration depth of 
the particles in a direction which connects the particle's 
centres of mass, pointing from i to We also intro- 

duce the penetration width of the particles defined as 



the direction perpendicular to the previous one. Thus 
we have 



•j(0 
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(6) 



valid for 



> 0. Let Cju-\ be a vector which defines 

a point Cjty of the application of the repulsive force, 
and which is taken as the mass centre of the overlapping 
region (J5J as shown in Fig. ^ Taking into consideration 
the fact the that particles have only spherical forms and 
collide when their overlap (j^J is positive we obtain 
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e Cj(0- ( ? ) 



Thus, at the beginning of a collision we have 



j'(i) x i 



r m + n 



(8) 



and one can find an explicit time t h -^ where the over- 
lap (0 is zero. Above notation allows us to analyse mul- 
tiparticle collisions where an individual particle i collides 
with neighbouring particles Therefore many over- 

laps (JSJ) indexed j(i) on the particle i may occur. This 
allows us to formulate the motion equations in the right 
form. 

Next we introduce the relative velocity of a particle i 
and a particle at point Cj^\ as 



u 



tin 



M f U j(i) 



rot 



(9) 



where u.TK and u r ?*s are linear and rotational relative 

J CO J CO 

velocities and Sijr^, Sj^^ are branch vectors connecting 
the mass centres of particles i and with the point 
Cj(i) of application of the repulsive force. Note that 
above values are defined for the global system of coor- 
dinates (x,y,z). To change this to the local system of 
coordinates we need to use the scalar product of the base 
vectors (@J per the vector of the relative velocity. Calcu- 
lating the branch vectors we obtain the following depen- 
dencies in the local system (£, ?y, () as 



where 
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(12) 



In the global system of coordinates the branch vectors 
become 



'T 



3 j(i),i 



'T 



(13) 



Using Eqs. and |j5Jl we translate the linear and rota- 
tional relative velocities in the global system of coordi- 
nates to the local one as 



U ; 



■ "' are the relative velocities operating 
in the normal direction to the contacting surfaces as 



u 



rot 



e • u' 



(14) 



where u'fi (i) , u <m 



shown on Fig. and u. 1 " 1 



~lin ~lin 
U i3{iV U r,3(i) 



~l' rot 



T:rot T.rot 



denote vectors of the relative velocities 

acting in the tangential direction (rotated by J to the 
normal). Additionally we use the same translation as 
presented by expression (|14H for calculations uj,- l = e • Wj, 
= e ' * n oro - er ^° obtain the angular velocities 
for particles i and in the local system of coordinates 

If a collision between a particle and a wall takes place, 
the overlap 10 is defined as 



>b 



xb 



and also we have 
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(15) 
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which is valid for 
,6 



>b 



> 0. In this case the point of 



application c^, s is defined by the following formula 
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where 
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xb 



(17) 



(18) 



becomes a normal unit vector which is perpendicular to 
the wall. When a particle-wall collision begins we ob- 
tain Cj(i\ = x,; + r i eb Qj(iy Moreover, one can find the 
explicit time t^.-. when the overlap i|15|) is zero. Expres- 
sions (|9l — (|14|) . defined for a particle-particle collision, 
may be redefined in simple way for a particle-wall colli- 
sion when and ijJjU) are zeros and unit vectors are 
also redefined as explained in previous considerations. 
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For example a component of the branch vector (|llfl is 
redefined for a particle-wall collision and takes the fol- 
lowing form 



(19) 



We neglect here any additional expressions necessary to 
describe the particle-wall collision. The reader can do 
this very easily in the same way as explained previously. 

Summarising our considerations, we introduced the 
above mathematical description which is necessary for 
the formulation of the motion equations and is also nec- 
essary for some forms of the repulsive force, acting for 
both particle-particle and particle-wall collisions. 



B. Motion equations 

The molecular dynamics method requires a discrete de- 
terministic approach in order to the model motion of an 
individual particle. Note that the particle may collide 
or lose contact with other particles. Therefore in mo- 
tion equations we need to add or reject some forms of 
the repulsive force and/or the attractive force in order to 
simulate the particle dynamics more realistically. In this 
paper we neglect the attractive force and we will con- 
centrate only on the repulsive force. Aganist this back- 
ground, let us describe the motion of an individual par- 
ticle by the following two sets of equations 



(20) 



suitable for particle motion without any collision, and 




J'(») 

Er,b coll 

E M co " 

Eiyrfc coll 

j'M>j(*)#* 



I 



EM, 



(21) 



which takes into account multiparticle collisions. The 
above sets of equations exist simultaneously over time 
and are dependent on the detection of a contact and the 
administration of the repulsive force-overlap path dur- 
ing the contact. In both Eqs. (|20|l and (|21|l F; denotes 
an arbitrary force which extorts the motion of a parti- 
cle, M/ is an arbitrary torque, P°°'j is a collisional force 
composed of the repulsive and friction forces and acts be- 
tween a pair of colliding particles, is also the colli- 
sional force operating on a particle- wall collision, M?°'^ 
and M^°" are collisional torques definitively operating 
on particle-particle and particle-wall collision. 

We need to define some of the criteria necessary for 
handling the above two sets of equations over the time of 
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the calculations. It was evidently shown in the previous 
subsection that for the beginning of a collision the over- 
lap given by expression © or l|15|l is zero. Thus we have 
the impact phase. However, some of the criteria for de- 
termining when the collision ends are unclear. Correctly 
predicting the separation time of two colliding particles 
is crucial in the calculation. Most papers assume the par- 
ticles separate at the time when the overlap returns to 
zero. As proved in |3^|, the repulsive force changes direc- 
tion at the time when the overlap returns to zero. This is 
contrary to experimental evidence, also shown and com- 
pared with some models in |35j. when the force does not 
change direction. An attractive force operating in oppo- 
site direction to the repulsive force has different origins 
and is not taken into account here. 

At this crucial point of our considerations, we need to 
introduce some definitions in order to predict correctly 
the beginning time of a particle collision and the time 
when the collision ends. Let us consider the time of cal- 
culations t G (0, T) where T represents the total time 
in which the calculations are performed. We also define 
the time step At in which we trace the system dynamics. 
Follow on from previous explanations we start with some 
conditions. 

Definition 1 //, within a time interval (t, t + At) detects 
the beginning of a collision between a pair of particles is 
detected then the overlap should fulfil the following 
conditions 



< and 



and therefore 



C m (t + At) 



> o 



0, 



(22) 



and then time tju\ G (i, t + At) is the time when the 
collision starts. 

Definition 2 If, within a time interval (i, t + At) the 
end of a collision is formed then the overlap J3|) and the 
normal component of the repulsive force R^j(i) should 
fulfil the following conditions 



> and 



C,m(* + Af) 



< 



d R ai;l) (t) > Q and Rcj(i)(t + At) > (23) 



and therefore 



£j'W (5(0 
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Ci W (t) 



> and 



Cj(i)(t + At) 



> 



id R Cj(l) (t)>0 and R Cj(l) (t + At)<0 
and therefore R^j(i) = 0; 



(24) 



and then time t^u-, G (t, t + At) is the time when the 
collision ends. 

In formulae (|23|l and l|24|) R^j(i) represents a normal com- 
ponent of the repulsive force. In the next subsection we 
will introduce a definition of this force. 



Definition 3 If, within a time interval (t, t + At) the 

overlap and the normal component of the repulsive force 
R^j(i) behave as follows 

t m (t)\\>0 arid ||c iW (t + Ai)|| >0 (2g) 
and R (j ( i) (t + At) -» 0+, 

then time tj,~ = At + t is the time when the collision 
ends. 

Definition 4 When the condition \2ty is fulfilled then 
linear and rotational components \1J$ of the relative ve- 
locity predict the following states after the collision: 

• rebound of particles without particle deformations 
f° r «<&<) (5 W ) ¥= and u l l n m (t e j({) ) has an op- 
posite direction (sign) to u^™^) (tj(i))> 

• torsion for w (i (t^J or sliding 

f° r K m feo) + or rollm 9 f° r K% (%)) * 

of particles without particle deformations for 

s&m (%)) = °' 

• non-permanent or permanent stick of particles 



without particle deformations for u 



t<) 



Stf w (%)) = and for u' t % (t^) + 

Definition 5 When the condition \24\j is fulfilled then 
components \1J$ of the relative velocity predict the same 
states as described by definition 4 but particle deforma- 
tions are noted. 

Definition 6 When the condition \25)) is fulfilled then 
the normal component 2^™^ y^jfj^j °f the relative veloc- 
ity predicts adhesion-induced plastic deformations of 
particles or breakage of particles depending on the hard- 
ness of the contacting surfaces. 

On the base of previous assumptions and definitions 1 
and 2 we introduce the collisional time between a pair of 



contacting particles as t c °^ — t e -^ 



r 



This time is 



determined by conditions (|23|l and l|24|l simultaneously 
In other words, when the overlap changes sign faster than 
the repulsive force changes direction or vice versa then 
the collision is finished. If particles are still in contact 
then the total contact time is significantly greater than 
the collisional time. If particles are separated then the 
total contact time equals the collisional time. As pre- 
sented in the first section, the collisional process com- 
poses the impact phase, the contact phase and the last 
phase formed after the contact phase. Moreover, when 
the formulation of the first and the last phases is infinites- 
imally short in time then the collisional time is predicted 
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by the contact phase. The contact phase is predicted by 
the repulsive force-overlap path. The adhesion or cohe- 
sion states extend the contact phase over time to infinity. 
In our approach adhesion and cohesion are eventually 
formed after the impact and they represent completely 
different phenomena which definitely result from the col- 
lision process. Generally, when we model impact dynam- 
ics we need to consider the balance between the repulsive 
force which is a direct reaction to the impact, and the 
attractive forces which are a result i.e. cohesion of par- 
ticles. Therefore, our collisional time t c °(t) a lso becomes 
the time of relaxation in which the collision process is 
stopped and novel states are formed. Most papers ne- 
glect this fact and identify the total contact time, which 
may increase to infinity, as the collisional one. 

Extending our considerations we notice that defini- 
tion 4 is suitable for the elastic collisions of particles 
because there are no deformations in contacting parti- 
cles - the overlap tends to zero faster than the repulsive 
force changes direction. In definition 5 we observe the 
opposite situation - the repulsive force changes direction 
faster than the overlap tends to zero. In our approach the 
collision process is fully controlled by the repulsive force 
except for the situation presented by formula (|25|l in def- 
inition 3. On the basis of definition 6, which results from 
definition 3, we are able to explain that the local stresses 
associated with deformations of contacting particles be- 
come sufficiently large so as to exceed the elastic limit 
of the materials to a result plastic flow occurs and 
the behaviour of particle adhesion differs from that pre- 
dicted by the elastic deformation theory |2l| ■ Adhesion- 
induced plastic deformations of contacting materials are 
evidently shown in some experiments 28] . Therefore we 
have two possible states resulting from the impact: par- 
ticle clusterisations when the colliding materials are soft, 
and fragmentation of particles when the colliding mate- 
rials are hard. 

Summarising this subsection we formulated two gen- 
eral forms of the motion equations and discussed precisely 
how to handle them. 



with sliding friction is for colliding particles which have 
different angular velocities in the normal direction and 
different linear velocities in the tangential plane. Sliding 
friction happens when slipping occurs in colliding parti- 
cles. When the relative linear velocity of the particles 
in the tangent direction reduces to zero, sliding friction 
is replaced by rolling friction. If the external forces are 
sufficiently small, the rolling friction reduces the veloc- 
ity until particle motion stops and static friction occurs. 
Considering the impact dynamics, we implemented the 
following mechanism in general form: torsion with slid- 
ing friction can change to rolling friction and the rolling 
friction tends to static friction. More details concerning 
the modelling of torsion, sliding and rolling friction can 
be found in j^, l3fij | . Here we show a description of the 
collisional force in global coordinates (x,y,z) as 
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(26) 



where P$") is the force acting in a static friction state, 
Pj°i) is the force occurring in a rolling state and 
is the force coupling the torsion-sliding state. The em- 
phasis in this paper is on the impact dynamics the static 
friction is only implemented in a simple form. A more de- 
tailed model of the static friction state requires analysis 
of the tangential displacement and possibly the inclusion 
of time dependent effects. According to Fig.^we need to 
define the collisional force in the local system of coordi- 
nates (£, rj, £)• Using a matrix of the base vectors Q) we 
introduce transition from the local system to the global 
ones as 

Tjsli _ T T}' sli r>rol T t>' rol 1 07\ 

where Pjf,-\) ^ are forces defined in the local system 
of coordinates as 



C. Collisional forces, collisional torques and the 
fractional interaction law 

With regard to motion equations l|21|) we introduce 
a mathematical description of the collisional forces and 
torques occurring in such a system. In the normal direc- 
tion to the contacting surfaces we apply only a repulsive 
force, completely neglecting any attractive forces. In j3f| 
one can find some forms of attractive forces and their 
physical meanings. In a tangential plane we introduce 
a system of friction forces and torques. According to the 
friction mechanism, the tangential friction force is one 
of four types: torsion with sliding friction, sliding fric- 
tion, rolling friction or static friction. Torsion friction 
occurs when colliding particles differ by their angular ve- 
locities in the normal direction w^i and oj^juy Torsion 
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In expression T sh 



rpsli 

nj(i)i 



rprol 



T 



rol 



represent 



components of the friction force in a plane (£,77) for 
torsion-sliding and rolling states, R(j(i) is a sum of the 
normal components of attractive and repulsive forces op- 
erating during a collision. As assumed in this paper, we 
neglect attractive forces and concentrate only on forms 
of the repulsive force. Some forms of the attractive forces 
can be found in 0, but the most well-known forms of 
the repulsive force are in P. Il3l Isflj . 

On the basis of preliminary results |15| we now intro- 
duce a model of the repulsive force in general form called 
the fractional interaction law. Thus we have 
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(29) 



where c^ya, are damping and spring coefficients with degree of impact energy into viscoelasticity of the mate 



the same meaning as in the linear interaction law rial and t T> a i {i) ( 

* ,1 l i n i i r i m 3'(«) 3'(«) V 



represents the overlap defined by formula JSJ, 

^(i) are start and stop times of a collision (not a to- 
tal contact) predicted by several definitions in the previ- 
ous subsection, as explained in 0] a.j(j\ is the conversion 



represents general form of the 

differential and integral operator of fractional order. Ac- 
cording to fractional calculus |U we introduce the 
definition of this operator in the following form 
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where t denotes actual time of calculations t G 



^ifil'^ifil)' ^ ne sum represents the initial conditions, 
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is the Caputo fractional derivative 
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where n^) = [a!j(i)] +1 and [•] denotes an integer part of Liouville fractional integral 
P ? w ( Cj(i)(t) ) is the Riemann- 



a real number, and t b T t 
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and PjU) — —oiju). Eqn. I|29|l represents the form of 
the repulsive force acting in the normal direction to the 
contacting surfaces. 

Now we introduce additional definitions of forces op- 
erating in the tangent plane. Here we define the normal 



force as N' ju) = [0, 0, R^jU)] ■ According to p| we define 
the friction force which is coupled between torsion-sliding 
friction as 



t ;« = -/* ( 
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where the friction coefficient is 

; lin 



ft 



M<i + (Ms - Md) e II tj(i) l 



where a is a numerical constant, fi s and /Zd are static and 
roo\ dynamic coefficients of friction. Moreover in formula H33|) 
the function T (Xj^)) is defined according to p| as 



(34) 



4 (A j 2 w +l)i;(A 3W ) + (A j 2 w -l)K(A 3( , ) ) 
3 TrXA(i) 



for A < 1 
for A > 1 



(35) 



where iC(Aj(j)) and E(Xj^) are the complete elliptic in- 
tegral functions of the first and the second kind, XjU) is 
the dimensionless quantity defined as 
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(36) 



The limiting values of the function T (Aj^)) are -^(0) = 
for torsion without sliding and lim J- (Aj(j)) = 1 for 

A j(i) ^oo 

sliding without torsion. 

According to [3(| we define the rolling friction force as 



n' vol 
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(37) 



where 



and 
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(38) 
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(39) 



The above expressions are necessary for the definitions 
of some collisional torques. Therefore we have the col- 
lisional torque operating from particle i to particle 



as 
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where M?'!,.^ is the coupled torsion-sliding torque, 



M 



rol 



represents the coupled torsion-rolling torque. 



Note that transition from the local system of coordinates 
to the global ones reads 



M ^w= e ^)-( M ^) + M i 

•m ■ (m; 



tor 

m 



-\j\rol 



rol 



tor 



(41) 



and according to Q we obtain 



irlx 



u 



tj(i) 



AT Ci(i) sign (w c< - oj (j{i) ) 



(42) 



We define the torsion torque as 



where the function T (^jU)) reads 



( 4 - 2A J 2 w) g ( A J w)+( A J 2 (.)- 1 ) g ( A J w) 



r ( A ^))- / 4 - 2A ^)K^7V( 2A ^-^ifc)^(^y) 



7rA, 



for Aj (i ) < 1 



for Aj(i) > 1 



(43) 



The limiting values of the function T (Aj/j)) are T(0) = 
| for torsion without sliding and lim TfAj-(i)) = 

A i(j) -»oc 

for sliding without torsion. Moreover, we introduce the 
sliding torque as 



tl/t' sli _ ' v T 1 ' sIi 



(44) 



Using an idea included in |3a ] we determine the rolling 
torque as 

M ^)^ s «(.)^w+ ss y(.) x %). ( 45 ) 



(|21|1 in order to simulate the dynamics of multiparticle 
collisions. The above system is mathematically complex, 
and therefore requires a numerical approach. An accu- 
rate solution to this problem was obtained by integrat- 
ing the system of ordinary differential equations (|20|l for 
particles moving individually by using Numerical Recipe 
routines |2fi| . Tracing the motion of individual particles 
over time we need to detect particle collisions in order 
to take into account collisional forces and torques in the 
system of differential equations. Using results presented 
in 0, 0, we have chosen the linked cell method to 
detect a collision. 



where ss i x is the torque created on the penetra- 
tion width JHJ. As noted in j3(|, the torque ss^ .^x x 
exists because the contact between two particles is not 
a single point but, due to deformation of both bodies, is 
a finite area. 

Summarising this subsection we determined a full de- 
scription of the forces and torques occurring in a collision. 
We neglect here a mathematical description of the colli- 
sional force P b S°J l and torque M hc ?'( acting between the 
particle- wall because one can easily produce these formu- 
lae taking into account i-jU) = 0, UjU) — 0, etc. in the 
above expressions. More details concerning particle-wall 
interaction can be found in [Tl| . 



III. SOLUTION PROCEDURE 

Throughout this section we will show how to handle 
the system of ordinary differential equations 1|20JI and 



Note that during particle collisions we need to solve 
the system H21(l where the fractional interaction law Ij29(l 
occurs. In this case we have a system of ordinary differ- 
ential equations with a mixture of operators: the integer 
derivative of maximal order equals two, the fractional in- 
tegral of order — <x,(i) and the fractional derivative of or- 
der ctj(j\ . Using fractional calculus [23|,|3(j we will present 
discrete forms of the fractional operators which are suit- 
able in our algorithm. Let us consider the duration of 
a collision over time t 6 (to,t nt ) where to represents the 
time when the collision starts and t n t is the time when 
the collision ends, nt denotes the division of the colli- 
sional time t into several time steps. Thus we obtain: 
h = tnt ~ t to , U = t + Ih, for I = 0, . . . , nt. If a function 
f(t) is constant within the step h then the discrete form 
of the Caputo fractional derivative becomes 
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T(n-a + 1) 



A x {t nt - t ) n a +J2( Al ~ " 



r 



i=2 



(46) 



where a £ 7?. + , n = [a] + 1 and [•] denotes an integer part 
of a real number, Ai = (£ z ) where is the deriva- 
tive of integer order n. Note that in formula (|46|l /(t) 



denotes the overlap JSJ). Taking the above assumptions 
into account we obtain the discrete form of the Riemann- 
Liouville fractional integral (|32[l as 



t iLf(t) = 



i 



r(/? + i) 



Si (t nt — to) 13 + ^ {Bi — -B;-i) (i«t — U-i) 13 

1=2 

I 



(47) 



where (3 S 7?. + and £>; = f (ti). The Discrete forms of 
the fractional operators makes it possibile to integrate 
the s yste m (|21[) by using any predictor-corrector proce- 
dure |26j with correction of the time step h. The correc- 
tion of the time step provides measures that allow us to 
determine the begin time when particles enter into a col- 
lision and the end time of particle collisions. It should 
be noted that the begin and end times are determined by 
several definitions presented in the previous section. 

Taking formulae l|26|) , (|40|l into account in calculations 
of particle contacts we need to find an accurate time 
needed to detect the switching between torsion-sliding, 
sliding and rolling processes. As described in the pa- 
per j3g, a simple way to calculate the switching time is 
to use a linear approximation method. 

Next we consider a problem occurring in the calcu- 
lations of friction forces l|33|> . (|37|) and the torsional 
torque l|42|) . When the relative velocity at the contact 
point changes from negative to positive or from positive 
to negative, it indicates that the signum function sign(a;) 
changes sign very fast in above expressions. This is not 
desirable as it influences the stability and convergence of 
the numerical calculations in a significant way. Therefore 
we modified the signum function introducing 



sign(x) = < 



£2-ei 



i . 

«2-ei ' 
1 



for 
for 
for 
for 
for 



x < -e-2 

-62 < X < -€l 
— El < X < €\ 

ei < x < e 2 
x > e 2 



(48) 

where x is the actual value registered during a contact 
(the relative velocity), t\, e 2 are numerical coefficients. 
This function is robust for x — > and gives a satisfactory 
result. 



IV. RESULTS AND THEIR ANALYSIS 

To illustrate the benefits of the fractional interaction 
law (|29|) in the dynamics of arbitrary multiparticle col- 
lisions, first we will demonstrate how this law operates 
in simple cases connected with a one dimensional prob- 
lem. First we will simulate a central collision between 
two particles. Fig.[2]shows the dynamics of a two-particle 
collision, which is represented by some variations in the 



overlap ||£|| ©, the linear relative velocity ( — 



and the repulsive force Rq (|29|l over time for different 
levels of the conversion degree a. Here we neglect the 
index j(i) because only two particles collide. Moreover, 
all vectors are converted to scalar values when a one di- 
mensional problem is considered. For the figure we have 
m eff = = 7 ' 06858 • 1(r6 fc 3> n = r 2 = 3 • 10- 3 m, 

k = 5000 -f, c = 0.1 kg /s. The initial relative velocity 

is set at £ = 0.5 — and three groups of variations in the 
conversion degree are taken into account. The first group 
is for a < (left column), the second is for < a < 1 
(middle column) and the third represents a > 1 (right 
column). Within the range < a < 1 we observe that 
collisional time t c increases when a is increased. It should 
be noted that the collisional time is registered when the 
repulsive force reaches zero, as presented in several 
definitions in the previous sections. Therefore, the over- 
lap IICH has some values at the time when a collision 
ends and deformations of the particle surfaces are noted. 
Analysing the behaviour of the relative velocity C over 
time we notice that this velocity changes direction for 
small values of a, which means thatparticle rebounds 
dominate. When a increases we can see that the rela- 
tive velocity tends to zero, which means that particles 
stick together. In other words, if a = 0, no viscous term 
in Eqn. (|29|l may occur and all the impact energy must 
be due to elasticity. In this case the overlap reaches zero 
at the same time as the repulsive force reaches zero. If 
a = 1, on the other hand, the impact energy is transfered 




FIG. 2: Behaviour of the overlap (top), relative velocity (middle) and force (bottom) over time for the fractional interaction 
law. 



through the viscous term. 

Extending our considerations for a > 1 we observe (right 
column on Fig. |2J) that the repulsive force is not gener- 
ated and tends to zero for t c — » oo, and therefore the 
overlap increases to high and unrealistic values. More- 
over, the relative velocity does not change direction and 
particles undergo the next time steps of the calculations. 
According to definition 6, presented in the previous sec- 
tion, the fragmentation of particles or permanent cohe- 



sion of particles is a direct result of the plastic flow of 
their contacting surfaces. The contacting surfaces are 
destroyed because deformations of contacting particles 
become sufficiently large so as to exceed the clastic limit 
of the materials, and we noticed particle clusterisations. 
This process is observed experimentally in |l2l |28( and 
may be modelled by the fractional interaction law i|29|) . 
Next we considered the behaviour of the overlap, relative 
velocity and repulsive force for a < (left column on 
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Fig- IHO ■ Larger negative values of the conversion degree 
a decreases the collisional time. The relative velocity 
changes direction but at the end time reaches larger ab- 
solute values in comparison to the initial relative velocity. 
As this is unrealistic all the solutions for a < are not 
taken into account. The aim of this example is to show 
the power of fractional calculus where more solutions are 
obtained in comparison to classical differential and inte- 
gral operators having integer order. However, we need 
to choose which solutions obtained by fractional calculus 
are suitable physically. 

In Fig. we constructed several mappings for the rel- 
ative velocity-overlap (left), force-overlap (middle) and 
force-relative velocity (right) where a changes from neg- 
ative to positive values. Analysing these mappings we 
found a set of criteria necessary to predict different states 
of particle collisions included in definitions in the previ- 
ous section. It should be noted that small positive values 
of a predict particle rebounds when particle deformations 
are practically neglected. When a tends to unity we also 
observe particle rebounds but particle deformations are 
visible and more energy is dissipated. As indicated in 
the left chart in Fig. when a is above unity the repul- 
sive force is not generated and this indicates instability in 
particle collisions. This instability takes the form of par- 
ticle fragmentation or permanent clusterisation of parti- 
cles after the collision. Therefore the conversion degree 
a is a ratio of the impact energy over the specific energy 
needed for the destruction of particle surfaces. This as- 
sumption should be validated experimentally, and this is 
the aim of our future investigations. Note that when the 
physical properties of colliding granular materials and the 
impact energy are fixed we still observe different values of 
energy dissipation after the collision. This can be easily 
seen when we compare the particle collisions for smooth 
surfaces of particles and for rough ones. The fractional 
interaction law can simulate this because the conversion 
degree a can change. 

In order to compare the fractional interaction law with 
other interaction laws, changes over time of the over- 
lap, the relative velocity and the repulsive force for two- 
particle collision were presented. We assumed parame- 
ters of colliding particles to be n = = 3 ■ 1CP 3 m, 
m eff = 7.06858- 1CT 6 kg, ( = 0.5 Moreover, we as- 
sumed the collision time between two colliding bodies the 
t c = 1CP 4 s and the restitution coefficient e r = 0.5. These 
assumptions are necessary to calculate the set of coeffi- 
cients required by different interaction laws, depending 
on the type of interaction law chosen. In Table [I] we list 
all the coefficients. Some of the expressions applied to 
calculate the coefficients for linear Q and hysteretic |3^| 
laws can be found in J25j. The formulae of the coeffi- 
cients used in the linear interaction law assumed that at 
the end time of a collision the overlap is zero. In Table [I] 
the "linearl" represents the above case. We assumed that 
the repulsive force reaches zero at the end time of a col- 
lision. Thus we have a set of coefficients called "linear2" 
also used for the linear interaction law. For the non- 



TABLE I: Coefficients for colliding particle surfaces being de- 
pendent on the interaction law used. 



law 



coefficients 



linearl 

linear2 

non-linear 

hysteretic 

fractional 



k n = 7316^, c„ = 0.0979^ 
= 5225 ^, c n = 0.0981 ^ 



1392000 - 



S 

c = 33.885 -H. 
fei = 3924 fc 2 = 15697^ 
k = 5225 c = 0.297 a = 0.3197 



linear [13j and fractional laws we performed a numerical 
test to find the values of coefficients which allow us to 
keep the assumed collision time and the restitution coeffi- 
cient in a two-particle contact. It should be noted that we 
obtained many sets of coefficients for the fractional inter- 
action law. Therefore for this law we establish the spring 
coefficient which has the same value as for the linear in- 
teraction law. Fig. 01 shows the behaviour of the overlap 
(top chart), the relative velocity (middle chart) and the 
repulsive force (bottom chart) over time where different 
interaction laws are taken into account. Analysing this 
figure we can confirmed that the interaction law fulfilled 
our assumptions concerning the collisional time and the 
restitution coefficient. It should be noted that the repul- 
sive force changes direction in the linear interaction law 
(linearl) for the set of coefficients calculated under the 
formulae found in |25| . This shows a deficiency in numer- 
ical calculations and should be rejected. Some changes 
in the values of the above coefficients give satisfactory re- 
sults in the linear interaction law (linear2). However, the 
repulsive force in the linear interaction law has a value 
at the beginning time which is independent on the set of 
coefficients used. This is also unrealistic behaviour in the 
linear law. 

Using different interaction laws we observed different 
overlaps at the end time of a collision. The highest value 
of overlap is for the hysteretic law and decreases for the 
fractional through the non-linear to the linear one. Note 
that for the fractional law we can find another set of co- 
efficients in order to fulfil our assumptions and to obtain 
another value of the overlap at the end time of collision. 
When we determine all parameters necessary to describe 
the dynamics of particle impacts then we obtain some 
values of the collisional time and the restitution coeffi- 
cient for this case. However when we still keep the above 
parameters and increase or decrease the surface rough- 
ness of colliding particles then we obtain values of the 
collisional time and the restitution coefficient differing in 
comparison to the previous values. As we did not change 
physical properties of this granular material therefore we 
have to save the steady value of the spring coefficient in 
all interaction laws. Changing only the damping coeffi- 
cient in the linear and non-linear laws and the unload- 
ing slope &2 in the hysteretic law we do not have any 
guarantee that we will obtain accurate values of the col- 




FIG. 3: Mapping the relative velocity-overlap (left), force-overlap (middle) and force-relative velocity (right) for the fractional 
interaction law. 
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FIG. 4: Comparison of the overlap (top), relative velocity 
(middle) and force (bottom) over time for different interaction 
laws. 



lisional time and the restitution coefficient reflecting the 
above cases. This is a disadvantage of the well-known 
interaction laws. In the fractional interaction law we 
have an additional parameter called the conversion de- 



gree a which causes some changes in the collisional time 
and the restitution coefficient. However, this requires 
some experimental data involving the impact dynamics 
of smooth and rough particles. These data will provide 
measures that allow some links to be made between the 
experiment and the coefficients of the fractional law. 

In order to verify the validity of the interaction laws 
for multiparticle collisions, the energies dissipated at each 
contact were compared. Here we introduce a measure of 
energy dissipation during multiparticle collisions which is 
the ratio of the kinetic energy evaluated in time over the 
initial kinetic energy. We define the total ratio of energy 
lost through multiparticle collisions as 



i=l 

nc 

E m i (i° 

i=l 



(49) 



where the superscript refers to the initial kinetic energy 
examined at time t = s and nc is the total number of 
colliding particles. 

We used a set of particles np vertically stacked over a bot- 
tom plate as shown in [la . l25j . We assumed the fol- 
lowing conditions ri = 0.0015 m, nrii — 1.41 • 10 -5 kg, 
ii = -0.5 for i — l,...,np. Gravity is set at zero. 
Taking into account the results presented by we cal- 
culated the energy dissipation as a function of the num- 
ber of considered particles np, which becomes the num- 
ber of colliding particles nc when at the begin time of 
the collision the distances between spheres equal zero 
I'j = Orn, for j = l,...,nc. Note that j — 1 repre- 
sents a collision between the first particle and the bot- 
tom plate and j = nc is a collision between the topmost 
particles. We also assume the collisional time between 
two colliding bodies t c — 10 -4 s and the restitution co- 
efficient e. r — 0.945. These assumptions are necessary 
to calculate some coefficients depending on the type of 
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interaction law chosen. The coefficients represent a colli- 
sion between two particles or between a particle and the 
bottom plate, where the plate mass is infinite. 
Fig. [S] shows the energy dissipation as being dependent 
on the number of collisions nc for different interaction 
laws used in the molecular dynamics method and also in 
the event driven method pi |22|. For linear, non-linear 




nc 

FIG. 5: Energy dissipation during multiparticle collisions for 
different interaction laws. 

and hysteretic laws we noted the same dependencies as 
in 0, [25| . This means that the "detachment" effect oc- 
curs. First, we considered the fractional interaction law 
for a steady value of the conversion degree ctj^ — 0.0258, 
for all binary collisions. In this case we obtained similar 
results for the hysteretic and fractional interaction laws. 
Thus the "detachment" effect also occurs in the fractional 
interaction law for the steady value of ctjny As written 
in [Tsl ] the kinetic energy obtained from the event driven 
technique is dissipated totally for nc ■ (1 — e r ) large. It 
should be noted that the basic interaction laws are valid 
for a two-particle collisions which are completely inde- 
pendent of other collisions. However, in multiparticle col- 
lisions we need to include mutual dependencies between 
several binary collisions. Taking this fact into account, 
we can obtain satisfactory results when the conversion 
degree cx-jU) changes in relation to the number of collid- 
ing particles. This was explained more precisely in |l5j . 
Therefore we propose a (nc) ~ 1 + exp(— nc) in order 
to keep a qualitative agreement with the event driven 
method. It should be noted that we cannot estimate cor- 
rectly a (nc) by direct comparison with the event driven 
technique. We require experimental data involving mul- 
tiparticle collisions. This data will provide measures that 
allow some links to be made between several coefficients 
in the fractional interaction law and the experiment. 

The last example simulates the dynamics of five par- 
ticles in three dimensional space for two values of the 
parameter a. The first value a = 0.01 indicates the 
strong repulsive state, i.e. particles rebound almost 
without dissipation of their energy. The second one 



for a — 0.97 represents the weak repulsive state where 
most of the impact energy is converted into material 
viscoelastcity. In the real behaviour of granular mate- 
rials we can easily observe such states, when we con- 
sider the collisions for contacting particles with smooth 
surfaces and for rough ones. For this simulation we 
assumed the following conditions r\ — 0.02 m, r 2 = 
0.01 to, r 3 = 0.007 to, r 4 = 0.005 to, r 5 = 0.009 to, 
6l = Q 4 = 2000 |£, g 2 = g 3 = g 5 = 1000 fa, xi - 
[0.0,0.1,0.23] m, x 2 = [0.001,0.125,0.205] m, x 3 = 
[-0.002,0.090,0.198] to, x 4 = [-0.004,0.120,0.186] to, 
X5 = [—0.001,0.1,0.18] to. Moreover, we consider a situ- 
ation where a particle with an initial linear velocity Ui = 
[0, 0, —5] y collides at different moments in time with 
particles which initially do not move (uj = [0, 0,0] — , for 
j = 2,..., 4). Particles do not rotate initially (u>i = -), 
gravity is set to zero and k = 1000 c = 1 for each 
pair of colliding particles. We also simplified values of 
the friction coeffcients putting into Eq. I|34fl a = and 
|U S = 0.5 for each pair of colliding particles. Fig. [S]shows 




FIG. 6: Behaviour of particle trajectories depending on strong 
(a = 0.01) and weak (a = 0.97) repulsions. 

the trajectories of the mass centres of five particles in 
three dimensional space for strong and weak repulsions 
as a reaction to the impact dynamics. The particles are 
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FIG. 7: Linear and angular velocities of particle "1" over time for strong (a = 0.01) and weak (a == 0.97) repulsions. 



numbered from "1" to "5". This simulation does not 
reflect the real motion of particles because we neglect ex- 
ternal forces, i.e. the gravitational force. We can only 
show how the fractional interaction law operates in the 
above conditions as being dependent on the conversion 
degree a. In the strong repulsive state (a = 0.01) we ob- 
serve linear particle trajectories. As a is increased and 
reaches the weak repulsive state (a = 0.97) we noticed 
different particle trajectories in comparison to the previ- 
ous state. According to the results presented in Fig. [21 
we can say that duration over time of the repulsive force, 
which is longer over time for higher values of a, has a sig- 
nificant influence on the particle trajectories. 
In order to explain more precisely what happens to par- 
ticle trajectories in strong and weak repulsive states, the 
velocities of one individual particle were analysed. Fig.0 
shows in global coordinates (x, y, z) the linear and angu- 
lar velocities of particle " 1" over time. In this figure the 
dashed lines represent particle velocities in the strong re- 
pulsive state, whereas continous lines indicate the weak 
repulsive state. We can observe clear jumps in particle 
velocities over time for the strong repulsive state. This is 



a result of the duration of a collision determined by the 
collisonal time between a pair of contacting particles. In 
this state we can notice binary collisions because several 
collisional times between the different pairs of contacting 
particles have shortest values in comparison to their sep- 
aration times, where particles move individually. How- 
ever, in the weak repulsive state we observe continous 
changes in particle velocities without the distinction of 
any jumps. This means that several collisional times be- 
tween the pairs of contacting particles overlap each other. 
The binary collisions are not distinguished here. 
Moreover, we analysed, in the local system of coordinates 
(Cj^iOj the angular velocities over time of particle "1", 
which collides with the particle " 5" . In the strong repul- 
sive state we observe smaller values of and u) v (these 
velocities are angular velocities predicted in the tangent 
plane as shown in Fig. ^) in comparison to the weak re- 
pulsive state. This means that torsion-sliding friction 
dominates in the strong repulsive state, where binary col- 
lisions are noted. In the weak repulsive state we observe 
that the angular velocities uj£ and oj v have higher val- 
ues than in the strong repulsive state. Thus we expect 
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the torsion-rolling state between particles " 1" and " 5" . 
However multiparticle collisions are noted in the weak 
repulsive state. 

In order to prove where binary or multiparticle colli- 
sions occur, some distributions of collisional times over 
the duration time of calculations are presented. Fig. |H1 
presents the sequence of segments of collisional times 
over the time of observation for strong and weak repulsive 
states. Continous segments represent collisional times for 
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FIG. 8: Sequence of collisional times depending on strong 
(a = 0.01) and weak (a — 0.97) repulsions. 

weak repulsion whereas strong repulsion is denoted by the 
dashed segments. Each segment representes one binary 
collision between a pair of contacting particles, i.e. "1"- 
"3" means the collision between particle "1" and particle 
" 3" . Analysing this figure we observe longer collisional 
times for the weak repulsive state in comparison to the 
collisonal times for the strong repulsive state. Moreover 
the collisional times ovelap in the weak repulsive state, 
therefore multiparticle collisions occur. 



dimensional space. We introduced a novel mathemati- 
cal description of this method which takes into account 
the division of the collision process into an impact phase, 
contact phase and another phase formed after the con- 
tact phase. We assumed that the impact phase and the 
phase formed after the contact phase are infinitesimally 
short in time. We redefined the collisional time so that it 
is predicted by the repulsive force-overlap path. On the 
base of preliminary results 0] we proposed an expres- 
sion for the repulsive force formulated under fractional 
calculus. The force can control the energy dissipation 
and the collisional time for an individual particle col- 
liding with many other particles. In multiparticle colli- 
sions we included the friction mechanism needed for the 
transition from coupled torsion-sliding friction through 
rolling friction to static friction. Therefore our model in- 
cludes multiparticle collisions in arbitrary forms. Using 
the fractional interaction law one can determine different 
states of particle repulsions, i.e. strong and weak repul- 
sive states. In the strong repulsive state binary collisions 
dominate, and torsion-sliding friction is the main friction 
mechanism. However, within the multiparticle collisions 
rolling friction is observed to be much stronger. 
The presented numerical results can be used to realisti- 
cally model the impact dynamics of an individual parti- 
cle in a group of colliding particles. In order to tune the 
model coefficients we require experimental data involving 
multiparticle collisions. This data provides measures that 
allow some links to be made between several coefficients 
in the fractional interaction law and the experiment. 
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